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Abstract 



Finding the densest sphere packing in <i-dimensional Euclidean space K is an outstanding funda- 
mental problem with relevance in many fields, including the ground states of molecular systems, col- 
loidal crystal structures, coding theory, discrete geometry, number theory, and biological systems. 

Numerically generating the densest sphere packings becomes very challenging in high dimensions 

-t— > 

S 

rJ-« ■ even for the restricted problem of finding the densest lattice sphere packings. In this paper, we 

Ch ■ 

apply the Torquato-Jiao packing algorithm, which is a method based on solving a sequence of 



due to an exponentially increasing number of possible sphere contacts and sphere configurations, 



linear programs, to robustly reproduce the densest known lattice sphere packings for dimensions 2 
through 19. We show that the TJ algorithm is appreciably more efficient at solving these problems 
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than previously published methods. Indeed, in some dimensions, the former procedure can be as 

o; 

If} , much as three orders of magnitude faster at finding the optimal solutions than earlier ones. We also 



study the suboptimal local density-maxima solutions (inherent structures or "extreme" lattices) to 
gain insight about the nature of the topography of the "density" landscape. 
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I. INTRODUCTION 



There has been great interest in understanding the packings of hard (i.e. nonoverlapping) 



particles because they serve as useful models for a variety of many-partic 



a variety oi i 
iquids |l|, U, 



e systems arising 
in the physical and biological systems, such as liquids [l|, |2j, glasses [3H5J], crystals jg48j, 
granular media |9l-ll2| . and living cells 13[. One outstanding problem is to find the densest 
packing of identical spheres in ci-dimensional Euclidean space M. d . This seemingly simple 
problem has proved to be a challenge for all but the most simple systems; it was not until 
2005 that a proof was successfully presented to confirm the centuries-old Kepler conjecture 
[14J ]. which states that the densest packing of spheres in three dimensions is the face-centered 
cubic lattice. For d > 4, there are no proofs for the densest sphere packings, although for 
d = 8 and d — 24 they are almost surely the Eg and Leech lattices, respectively 15]. 
Interestingly, these two lattices have also been used to construct 10- and 26-dimensional 
string theories, respectively jig, Il7|. 

In recent years, high- dimensional dense sphere packings have attracted the attention 
of physicists because of the insights they offer about condensed-phase systems in lower 
dimensions J5|, |l2J, ll8H20] . It is noteworthy that the general problem of finding the densest 



sphere packings in IR d (and other spaces) is directly re 



over communication channels resistant to noise 
geometry and number theory 



21 



22 



evant to making data transmission 



22] and of intense interest in discrete 



23j. The densest sphere packing problem is also deeply 



linked to the covering, quantizer, number variance, and kissing number problems, with which 



it shares the best known solutions in a variety of dimensions 



22 



24 



251 ]. Clever analytical 



methods have been used to discover dense packings in high dimensions (i.e., d > 4) but this 



approach becomes less efficient as d increases, especially because lessons 
dimensions cannot be used to construct packings in higher dimensions 22 



earned in lower 



26]. 



Numerical methods have only recently emerged to discover the densest pac king s in high- 
dimensional spaces. One such method devised by Kallus, Elser, and Gravel 19], is based 
on the "divide and concur" framework in which a dense arrangement of overlapping spheres 
is gradually relaxed until none of the spheres overlap. Another method formulated by 
Andreanov and Scardicchio 20] takes advantage of the fact that all densest lattice packings 



are also perfect lattices (defined precisely in Sec. IIVI) . which are finite in number 22]. The 
densest lattice packings can therefore be obtained by randomly exploring the space of perfect 



lattices. The efficiency of both algorithms plummets as d grows larger, preventing them from 
being effectively used in very high dimensions [27] . 



In the past twenty years, the Lubachevsky-Stillinger (LS) algorithm 



has served as 



a standard 
dimensions 



or generating dense packing of various shaped hard particles in two and three 



29 



31] . However, since the LS algorithm is based on a particle-growth molecular 



dynamics simulation, it is extremely computationally costly to use it to generate jammed 
dense packings with high numerical accuracy, especially as d grows beyond three dimen- 
sions. A recent improvement on the LS algorithm is the Torquato-Jiao (TJ) algorithm [321 ]. 
which replaces the molecular dynamics with an optimization problem that is solved using 
sequential linear programming. In particular, the density of a sphere packing (fraction 
of space covered by the spheres) within an adaptive fundamental cell subject to periodic 
boundary conditions is maximized. The design variables are the sphere positions (subject 
to nonoverlap), and the shape and size of the fundamental cell. The linear programming 
solution of this optimization problem becomes exact as the packing approaches the jamming 
point 12j. The TJ algorithm has been found to be a very powerful packing protocol to 
generate both maximally-dense packings (global maxima) and disordered jammed packings 
(local maxima) with a large number of identical spheres (per fundamental cell) across space 
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34]. 



dimensions [12J as well as maximally dense binary sphere packings 

In this paper, we specialize the TJ algorithm to the restricted problem of finding the 
densest lattice sphere packings in high dimensions. In a lattice packing, there is only one 



sphere per fundamental cell 35]. Even this limited prob 



em for d > 4 brings considerable 



challenges; its solution has been proven only for d < 8 49] and d — 24 1 5J], a nd it is closely 
related to the shortest- vector problem, which is of NP-hard complexity |36| . Additionally, 
most of the densest known sphere packings for d < 48 are lattice packings 22j, |26( . Tackling 
the lattice problem is thus a necessary first step prior to attempting to solve the much more 
complicated general problem of finding the densest periodic packings. A periodic packing of 
congruent particles is obtained by placing a fixed configuration of N particles where N > 1 
with in one fundamental cell of a lattice, which is then periodically replicated without 
overlaps. 

The outline of the rest of the paper is as follows: Sec. [Ill describes the implementation of 
the TJ algorithm for the special case of lattice sphere packings. In Sec. IIHI we motivate the 
choices that we make for the initial conditions and relevant parameters in order the various 



problems across dimensions. In Sec. IIVt we apply the TJ algorithm for 2 < d < 19, and show 
that it is able to rapidly and reliably discover the densest known lattice packings without a 
priori knowledge of their existence. The TJ algorithm is found to be appreciably faster than 



previously published algorithms 19j,|20|. We also demonstrate that the suboptimal-lattice 
solutions (i.e., the local maxima "inherent structures") are particularly interesting because 
they reveal features of the "density" landscape. In Sec. [V] we close with some concluding 
remarks and a discussion about possible improvements and other applications of the TJ 
algorithm. 

II. APPLICATION OF THE TJ ALGORITHM TO FINDING THE DENSEST 
LATTICE SPHERE PACKINGS 



The basic principle behind the TJ algorithm [32| resides in the fact that finding the 
densest sphere packing can be posed as an optimization problem with a large number of 
nonlinear constraints (such as nonoverlap conditions between pairs of particles) which can 
be solved by solving a series of linear approximations of the original problem. Its solution 
eventually converges toward a local or global optimum. While global optimality cannot 
be guaranteed, it has been shown that the TJ algorithm frequently reaches the globally 



densest packings [32J. The TJ algorithm was formulated for the general problem of finding 
dense periodic sphere packings. Here we describe its implementation for the special case of 
determining the densest lattice sphere packings, which reduces the problem to optimizing 
the shape and size of the fundamental cell, since no sphere translations are involved. It is 
interesting to note that the TJ algorithm can be viewed as a hard-core analog of a gradient 
descent in the space of lattices for energy minimizations for systems of particles interacting 



with soft potentials as described by Cohn, Kumar, and Schurmann [2 

Before explaining the numerical details of the TJ algorithm, we need to define some 
mathematical quantities. A (^-dimensional lattice A is composed of all vectors that are 
integer linear combinations of a set of d basis vectors mi, ..., m^, 

P = rii mi + n 2 m 2 H h n d m d , (1) 

where rij are the integers (j = 1,2, ... ,d) and we denote by n the corresponding column 
vector with such components. Using the generator matrix Ma, whose columns are the basis 



vectors, allows us to explicitly write the lattice set: 

A = {M A n : n G Z d } . (2) 

One useful property of Ma is that its determinant is equal (up to a sign) to the volume of 
the lattice fundamental cell. We can then write the lattice packing density as the ratio of 
the volume occupied by spheres of diameter D to the volume of the fundamental cell: 

m = «£&., (3) 

rV ; |detM A | 
where 

"« = WTW) (4) 

is the (i-dimensional volume of a sphere of radius R and T(n) is the Euler gamma function. 

The problem of finding the densest lattice packing of spheres in d dimensions can be 
expressed as: Find the d x d generator matrix Ma with minimal determinant | det Ma | , 
under the constraint that all non-zero lattice vectors MAn, neZ^ {0}, are at least as long 
as D. 

For this problem, the Torquato-Jiao algorithm consists of the following four steps: 

1. Randomly create a generator matrix Ma according to some stochastic process. 

2. For a given influence sphere radius Ri > D, find all of the non-zero lattice vectors it 
contains, i.e., compute {v = M_^n : n 6 Z d \ {0} A |v| < Ri}. 

3. Solve a linearized version of a problem, for which the objective is to maximize <fi 
(equivalent to minimizing | det Ma |) and the constraints are that none of the vectors 
calculated in step [2] become shorter than D. 

4. Consider whether the algorithm has converged to a lattice that is a stable maximum 
in (either the densest lattice packing or a local maximum inherent structure [371] ). 
If it is the former, repeat the procedure starting from step [2j If it is the latter, the 
solution has converged to a local or global optimum and the procedure is terminated. 

In what follows, we provide a more detailed explanation of these four steps. 

A. Initialization 

There are many possible methods to initialize the generator matrix Ma- Any candidate 
procedure must both satisfy the minimal length constraint and adequately sample the space 



of all lattices. The former is trivially satisfied by rescaling the matrix if the minimal length 
constraint is violated. In order to satisfy the latter condition, we mainly use Gaussian initial 
lattices, in which each coefficient of their generator matrix Ma is an independent normal 
variable N(0, a 2 ) with a variance a 2 . These matrices have the property that each of their 
lattice vectors (columns of Ma) have independent orientations with no given preference for 
any particular direction. To compare this against a different initialization method, we also 
consider initial lattices for which Ma is the sum of the generator matrix of a specific lattice 
packing (such as the d- dimensional checkerboard lattice Da or the hypercubic lattice Zj, see 
Appendix [X] for the definitions of these lattices) and one of a Gaussian initial lattice. 

B. Finding short vectors 

Finding all of the vectors for an arbitrary lattice that are within a small given radius Ri 
from the origin is a complex problem in high dimensions. Indeed, the problem of finding 
the shortest lattice vector for a given lattice A grows superexponentially with d and is in 
the class of NP-hard (nondeterministic polynomial-time hard) problems 



. One efficient 



39l . The influence sphere radius Rj can be 



method to solve this problem can be found in Ref. 
any value larger than the sphere diameter D, and may vary from one iteration to the next. 
It is found that the algorithm is largely insensitive to the value chosen for Ri, which is to 
be contrasted to the results for periodic packings, where larger Rj values favor the densest 



packings over inherent structures 32] . Since the computational cost of this and the following 



steps quickly increases with Ri, we opt to use the nearly minimal value Rj = LIZ?. 

C. Solving the linearized problem 

The only linearized problem variables in the case of the implementation of the TJ algo- 
rithm in the case of a lattice packing are the coefficients of the dx d symmetric strain tensor 
£ 40|. The modified generator matrix is then 

M A -)• M A + eM A . (5) 



The constraint that a vector originally at position v = Maii remains at least as large as D 
can then be written as 

n T M^M A n + 2n T M{eM A n + n T MX£ T £M A n > D 2 , 

v T v + 2v T ev + v T e T ev > D 2 . (6) 

This constraint is linearized by dropping the term that is quadratic in e: 

2v T ev > D 2 - v T v. (7) 

It should be noted that the term (v T £ T £v) that has been dropped is non-negative, which 
means that every set of variables that satisfies inequality (J7|) also satisfies inequality ([6]). This 
is different from the equivalent constraints for periodic packings, for which the quadratic 
term may be negative due to the interaction between the lattice deformation and the particle 
displacements. This avoids the necessity of either adding a constant term to the constraint 
or rescaling the system if spheres are found to overlap, which is the case for the general 



periodic packing problem [32 ] . 



Additionally, extra constraints must be added to prevent vectors that could be outside 
the influence sphere from becoming shorter than D: 

2v T sv > D 2 -Rj 
v T sv D 2 /R 2 - 1 _ 
^v" " 2 = " A ' (8) 

where the length of the vector has been chosen as its smallest possible value (Ri). A 
simple yet robust method to ensure that inequality (jSJ) is satisfied for all vectors outside of 
the influence sphere is to bound the lowest eigenvalue of e from below by —A. There are 
multiple ways to write linear constraints on e such that its eigenvalues are all larger than 
—A. One such way is given by 

< Diagonal element of e < oo, (9) 

— — — — — - < Off-diagonal element of e < — — — — . (10) 

Finally, the determinant of the modified generator matrix (assuming that det M A > 0) is 

det M A det (1 + e) = det M A (l + tr s + 0(e 2 )) , (11) 



where I is the d- dimensional identity matrix. The linearized density <j) is thus 

[1-tre], (12) 



n*J 



where 0o is the density for the initial generator matrix Ma and we used the fact that the 
density is inversely proportional to the fundamental cell volume. We can see from the above 
relation that maximizing the lattice density is equivalent to minimizing the trace of the 
strain tensor e. Unlike the linearized constraints ((71), ([9]) and ( TTOl) . which are conservative 
in that as long as they are satisfied the nonlinearized constraints will always be satisfied, 
the objective function (fl2"l) may have the wrong sign due to the nonlinear term having an 
unknown sign. In the situation where the updated lattice has a larger determinant than 
the original matrix, we halve e (multiple times if necessary) to ensure a lower updated 
determinant. This prevents the algorithm from oscillating between multiple lattices and 
forces it to eventually converge. 

D. Convergence criterion 

The algorithm is considered to have converged if the sum of the squared coefficients of 
e is below a small threshold value (10~ 12 for this paper). This is numerically equivalent 
to saying that all lattices in the neighborhood of the current lattice are less dense. This 
resulting lattice is therefore a local density maximum ("inherent structure" or "extreme" 
lattice, as elaborated in Sec. IIVBI) . Such a lattice is also strictly jammed, since any possible 
deformation requires an increase in the volume of its fundamental cell [III, |4jJ, |42 1 . 



III. STUDY OF PARAMETERS AND INITIAL CONDITIONS 

The ability of TJ algorithm to discover the densest lattice packings can potentially be 
affected by the influence sphere radius Ri, the lowest eigenvalue of the strain matrix A, and 
by the choice of the initial lattice. This section is dedicated to the study of their impact on 
the algorithm and to explain our choices for them in the following sections. 



The TJ algorithm is deterministic 43J, and therefore the initial lattice fully controls 
the resulting final lattice for given parameters Ri and A. For example, employing initial 
lattices that are very close to the known densest lattice, not surprisingly, results in a very 
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high success rate in obtaining that lattice. On the flip side, it would almost certainly never 
be able to discover a hypothetical denser lattice. It would therefore be misguided to use 
configurations that are near the known densest lattice as the initial conditions. However, 
allowing initial lattices that are very bad packers could result in a low success rate or a 
large convergence time for success. Thus, good choices for initial lattices involve a delicate 
balance between their diversity and an ability to relax quickly to dense lattices. 

TABL 



22, 



I: Frequency at which the densest known lattice packing in 13 dimensions, the K13 lattice 
26j, is obtained for various parameters using the TJ algorithm. For all sets of influence sphere 
radii and initial conditions, 10000 lattice packings have been generated, excepted for Rj = 2.0 
where only 3000 packings were generated. The calculations were performed on a single thread on a 
2.40 GHz processor using the Gurobi linear programming library ||44|. Since the run time strongly 
depends on the computer running the program and how well the code is optimized, it should only 



be used as a rough indication of the 


program efficiency. 








Sphere of influence radius 


Initial conditions 


Success rate (%) 


Average time per trial 


(sec) 


R r = l.LD 


Gaussian 


8.61 


5.0 




R r = l.LD 


Z?i3 + noise 


8.21 


5.5 




Ri = LID 


Z 13 + noise 


8.58 


5.2 




Rl = l.LD 


Invariant distribution 


8.08 


29.2 




Ri = 1.02D 


Gaussian 


8.53 


12.0 




flj = 1.5L> 


Gaussian 


7.61 


69.9 




R r = 2.0L> 


Gaussian 


6.87 


1938.5 




variable Rj, ~ 200 constraints 


Gaussian 


7.97 


6.3 




variable Rj, ~ 2000 constraints 


Gaussian 


7.95 


17.6 




variable Rj, ~ 2000 constraints, reduced A 


Gaussian 


8.58 


108.7 





Table [T] shows numerical results in 13 dimensions. The initial lattices are taken from 
four different distributions, using six different influence sphere radii. The TJ algorithm 
typically succeeds at generating the densest known lattice packing with a high probability. 
However, it has a relatively lower success rates for the cases d — 13 and d > 17. We thus 
purposely choose the 13- dimensional case to probe the best choices for the initial conditions 
and algorithmic parameters because of its abnormally low success rate in comparison to cases 



d < 16. Its low success rate results in better sensitivity to algorithm parameters compared 
with dimensions that have naturally higher success rates. Similar parameter dependence 
has been observed for other dimensions. 

The Gaussian initial condition, as previously explained in Sec. Ill A[ selects each coef- 
ficient of Ma from independent normal distributions with variances a 2 = D 2 . The initial 
conditions referred to as Dd + noise and Z d + noise starts with the generator matrices for 
the checkerboard Dd and hypercubic Z d lattices (these lattices are defined in Appendix [AT) , 
respectively, with nearest-neighbor distance equal to D plus some noise. Specifically, we 
add normal noise to each coefficient of Ma with a variance a 2 = D 2 /100. The final initial 
condition type that we attempt to employ, which we call an invariant distribution, generates 
the lattice from an approximation of the invariant lattice distribution, using the algorithm 
described in Ref. |45| with p = 10007. For all of these initial conditions, the nearest neighbor 
distance is calculated and the lattice is rescaled to avoid any sphere overlap. 

As can be seen in Table [U the different initial conditions that we have used result in 
similar success rates. We therefore use the Gaussian initial condition to generate the initial 
lattices for all subsequent calculations, since it lacks both the potential bias that the Dd 
+ noise and Z d + noise initial conditions share, and it does converge much faster than the 
invariant distribution. 

The main parameter influencing the efficiency of the TJ algorithm is the influence sphere 
radius Ri, which can either be fixed or vary from one iteration to the next. A radius that 
is too large leads to a large number of extra constraints for the linear program, greatly 
increasing its complexity. By contrast, if Ri is too close to D, then the constraints on the 
shear matrix s will be too restrictive [see Eqs. (jBJ), © and (TTOl) ]. This, in turn, only allows 
the lattice to deform very slowly, thereby requiring many iterations before convergence. A 
compromise between both is to use a variable Rj, such that the number of vectors inside 
the sphere of influence stays relatively constant, thus initially allowing a fast convergence 
when <fi is small, without needing numerous constraints when <fr gets close to its maximum. 
We use the following rough approximation to select Rf. 

Number of constraints ~ - — — , - , , (13) 

2|detM A |' v ; 

where the factor of one-half comes from the observation that for every vector v in a lattice, 
there is another one of identical length — v which does not need to be explicitly constrained. 

10 



A final parameter that can be modified is how much the lattice is allowed to deform at every 
iteration. As a test case, we divide the value of A by 10 to check whether an increased value 
of Rj provides benefits other than allowing larger strain matrices. 

From Table HI we can see that increasing Rj does not increase the success rate (it actually 
negatively affects it), while it significantly increases the run time. Therefore, the following 
calculations will be done using a small influence sphere radius of Rj = 1.1D. We attempted 
to adjust Ri as a function of dimension d to improve success rates for large d, but this 
proved to be fruitless. The radius Ri only weakly impacts the success rate, but its value has 
a dramatic influence on the time per trial, which gets multiplied by 400 when Rj is increased 
from 1.1D to 2.0 D. Therefore, one should decide on a choice of Rj so as to prioritize a faster 
execution speed over an increased probability of reaching the densest lattice packing. 

IV. RESULTS 

Here we describe the results we obtain by applying the TJ algorithm to find the densest 
lattice packings in dimensions 2 through 19. We compare our results with those obtained in 



previous investigations 19|, |20j . We also provide the frequency of time that the T J algorithm 



finds local versus the densest known global maxima. 

A. Finding the densest lattice packings 

We have applied the TJ algorithm for dimensions d = 2 through d = 19, and found 
the densest currently known lattice packing for each of them. The algorithm is robust in 
that it converges rapidly to the optimal solutions in most dimensions. Not surprisingly, 
except for the trivial d = 2 and d — 3 cases, it does not reach the optimal solution for all 
initial conditions. Therefore, even though the probabilities of finding the densest packing 
on the first attempt was high (greater than 19% for d < 12 and 14 < d < 16), we typically 
needed multiple trials (i.e., different random initial conditions) to guarantee that the densest 
lattice packings were among these. Consequently, the quality of such a global optimization 
algorithm is preferably measured using the time required per successful trial instead of simply 
the time per trial or the success rate. Table HT1 describes the rate at which the TJ algorithm 
produced the densest known lattice packings for dimensions d = 2 through d = 19 and the 
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TABLE II: Frequency at which the densest known lattice packing is obtained using the TJ algorithm 
for d = 2 through d = 19 together with the lattices packing fraction (j) and kissing number Z . The 
number of lattice packings generated is 10000 for d < 18 and 100000 for d = 19. The influence 
sphere radius Ri = 1.1D and the initial lattices are generated using the Gaussian initial condition. 
See Appendix [A] for the definitions of the various lattices. The comments in Table U concerning 
computational times also apply here. 



d 


Densest 
lattice packing 





Z 


Success rate 


Time per trial 
(sec) 


Time per successful 
trial (sec) 


2 


A 2 


0.9069 


6 


100 


1.7 x 10" 5 


1.7 x 10~ 5 


3 


D 3 


0.7405 


12 


100 


8.0 x 10~ 5 


8.0 x 10~ 5 


4 


D 4 


0.6169 


24 


74.31 


5.6 x 10~ 4 


7.5 x 10~ 4 


5 


D 5 


0.4653 


40 


97.41 


8.0 x 10~ 3 


8.2 x 10~ 3 


6 


Eq 


0.3729 


72 


89.72 


0.019 


0.022 


7 


Ej 


0.2953 


126 


91.91 


0.046 


0.050 


8 


Eg 


0.2537 


240 


84.16 


0.33 


0.40 


9 


A 9 


0.1458 


272 


43.82 


0.21 


0.49 


10 


Aio 


0.09202 


336 


22.74 


0.49 


2.1 


11 


K n 


0.06043 


432 


19.39 


1.1 


5.7 


12 


K l2 


0.04945 


756 


33.30 


2.7 


8.2 


13 


K 13 


0.02921 


918 


8.61 


5.0 


58 


14 


Au 


0.02162 


1422 


20.69 


10 


51 


15 


A15 


0.01686 


2340 


23.78 


16 


65 


16 


Al6 


0.01471 


4320 


22.50 


51 


227 


17 


A17 


0.008811 


5346 


1.65 


55 


3.4 x 10 3 


18 


Al8 


0.005928 


7398 


0.10 


79 


7.9 x 10 4 


19 


A19 


0.004121 


10668 


0.009 


162 


1.8 x 10 6 



average time required per successful trial. We determine whether we achieved the densest 



known packings primarily by comparing the packing density and the kissing number Z 
number of spheres that are in contact with any given sphere) with published data 22 



(the 



26] 



Additionally, we calculate theta series (the generating functions for the number of vectors 



12 



with specific lengths in the lattices 22J) up through the first few coordination shells. 



The time required by the TJ algorithm to generate the densest known lattice packings 



is appreciably smaller than the times reported in Ref. [19|: approximately 4000 and 25000 
seconds per successful packing for d — 13 and d = 14, respectively. The times required by 
the TJ algorithm of 58 and 51 seconds are orders of magnitude lower, indicating a genuine 
algorithmic improvement that cannot be attributed to the type of computer employed nor 
to implementation details. 

The authors in Ref. |20| do not state precise run times for all dimensions, but report 
that, after generating more than 10 5 lattices, their algorithm is unable to discover the 
densest known lattices for d = 14 through d = 19. Since generating 10 5 lattices using their 
algorithm takes at least several hours, the TJ algorithm's ability to successfully generate 
the densest lattice packings in minutes for d < 16 is a tremendous speed-up improvement. 
Using more computing power, the authors in Ref . |20| are able to reliably obtain the densest 
known lattice for d < 17 using their algorithm [46]. For example, their calculations took 
four days (~3x 10 5 seconds) for d = 14, which is three to four orders of magnitude longer 
than our own calculations (see Table HT1) . 

The fact that the TJ algorithm was unable to find any denser lattice packings than the 
densest known lattice packings reinforces the evidence that these are indeed the densest 
lattice packings for d = 2 through d = 19. Although this evidence is not as strong for d — 18 
and d = 19, due to the rare occurrences of the densest lattice packings, the evidence is quite 
strong for d < 17. 

One particular aspect of the success rates shown in Table [III is that they do not decrease 
monotonically with increasing dimension. Dimensions that are notably difficult are d = 4 
and d = 13, and neither case can be explained by lattice packings with unusual properties, 
since d = 5 and d = 12, respectively, share similar packings, but not the relatively low 
success rates. We will attempt to explain this phenomenon, along with the sharp decrease 
in success rates at d — 17, in the following section. 

B. Inherent structures 

The TJ algorithm is intrinsically a local density maximization algorithm. As such, it can, 
and often does, converge locally to the densest lattice packing associated with a given initial 
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configuration, i.e., an inherent structure 32], that are not necessarily the global maxima. 



These local maxima are analogous to the inherent structures of a continuous potential. The 
study of these inherent structures are of fundamental interest in their own right because they 
offer insight about the nature of topography of the "density" landscape and understanding 
the frequency of their occurrence could potentially lead to improvements on the algorithm. 

One interesting property of the density landscape associated with the lattice packing 
problem is that all of its inherent structures are extreme lattices, i.e., they are both perfect 
and eutactic |47J . Only a finite number of distinct extreme lattices exists for any dimension, 
which explains how the TJ algorithm is able to always reach the ground state for d = 2 and 
d = 3, for each of which only a single extreme lattice exists. However, as d increases, the 
number of extreme lattices grows quickly, possibly exponentially fast. It is thus remarkable 
that the TJ algorithm can reliably yield the densest lattice packing from the large set of 
possible end states. This indicates that the "basin of attraction" of the ground state is much 
larger than the basins of attraction of the local-maxima inherent structures. The relatively 
lower success rates for some dimensions (d — 4, d — 11, d — 13, and d > 17) can then be 
understood as being due to smaller than usual basins for the corresponding ground states. 
The cause of this reduction and whether the symmetry of the inherent structure is lower 
than that of the ground state or some other effect is still unknown and warrants further 
investigation. 

As seen in Table IHIl some inherent structures are degenerate in the sense that multiple 
lattices share the same packing density. A peculiar property that these degeneracies share 
is that their appearance rate is far from constant. For example, it goes from 9.24% for the 
A™ 11 to a mere 0.03% for the A^ ax . Since both of these are laminated lattices, why does one 
occurs more frequently than the other? One possible reason is that for all these degeneracies 
but one, the lattices with smaller kissing number are more likely to be generated. In the case 
of A^ m and A™ ax , their kissing numbers are respectively 624 and 648. This is consistent 
with previous work which has shown that for packings with many particles per fundamental 
cell, the TJ algorithm has a propensity to generate isostatic packings from random initial 
conditions, where the number of interparticle contacts is equal to the number of degrees of 
freedom of the problem [32J . 

Figure [1] shows that as the dimensionality increases, the inherent-structure densities tend 
to become concentrated around a specific value instead of being spread over a range of pos- 
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TABLE III: Second and third highest-density inherent structures (locally densest lattice packings), 
including their packing density (ft, kissing number Z, and success rate from the TJ algorithm. See 
Table [TT] to compare to the densest lattice packings. The number of lattice packings generated for 
each dimension is 10000 for d < 18 and 100000 for d = 19. Multiple lattices with equal density are 
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grouped together and written in ascending kissing number order. See Ref. 
of the following lattices: A§ , E£, P7.3, P7.5, Kg, Dimll (named dimllkis422 in the reference), 
K\ A , K\ A , h\ 5 , K\ 5 , k\ G , and K\§ . Lattices that were not identified in Ref. |26j and found here are 
denoted as U^, where n is used to distinguish different lattices at some fixed dimension d. 
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FIG. 1: Probability density functions for the packing density cf> (left) and probabilities for the 
kissing number z (right) of the lattice resulting from the TJ algorithm for (a) d = 13, |(b)| d = 15, 
and (c) d = 17. The minimal value of the kissing number Z m i n = d{d + 1) is 182 for d = 13, 240 
for d = 15, and 106 for d = 17. 



sible densities. This concentration tendency is caused by the rapid increase in the number of 
such low-density inherent structures for large d, which eventually overwhelms the algorith- 
mic bias toward high-density lattices. This explains the dramatic reduction in success rates 
in Table [III for d > 17. The kissing number has a similar behavior to the packing density, 
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resulting in the fact that most of the generated lattices for d > 17 have an identical low 
kissing number. Since these are locally-optimal solutions, a local deformation of the lattice 
would either decrease its packing fraction or makes the central sphere and its neighboring 
spheres overlap. Therefore, we can define a lower bound on the kissing number by exploiting 
the fact that, for a linear program to have a unique feasible solution, it requires at least 
one more active inequality constraint than the number of degrees of freedom. Since the 
problem possesses d(d + l)/2 degrees of freedom (the number of independent components 
of e), 1 + d(d + l)/2 active inequality constraints are required for the problem to be fully 
constrained. One of these constraints comes from the density being at a local maximum, 
while each pair of kissing spheres adds a single constraint. Consequently, the minimum kiss- 
ing number of a lattice inherent-structure in d dimensions is Z min = d(d + 1). Referring to 
Fig. [H we observe that as d increases, the proportion of generated configurations with a kiss- 
ing number equal to Z min increases rapidly relative to all other kissing numbers. Since the 
best known lattice packings have high kissing numbers (nearly the same or equal to highest 
known kissing numbers), the tendency of the TJ algorithm tendency to favor lattices with 
minimal kissing numbers further explains its low success rates for d > 17. 

V. CONCLUSIONS AND DISCUSSION 

In this paper, we have shown that the Torquato-Jiao algorithm is able to quickly find 
the densest known lattice packings for d < 19. The TJ algorithm is found to be orders of 
magnitude faster than the previous state-of-the-art lattice packing methods [l9j, |20j . This 
makes the TJ algorithm the fastest current numerical method to generate the densest lattice 
packings in high dimensions. 

While we limited our present study to d < 19, the TJ algorithm can be employed to gen- 
erate dense lattice packings in higher dimensions at greater computational cost. We expect 
that dimensions d = 20 and d — 21 would be manageable with more computing resources, 
but improvements to the algorithm would be required to study d > 22. One possible ap- 
proach to increase the likelihood of generating a dense lattice packing for d > 22 would be 
to include ad hoc methods in between the TJ-algorithm steps that favor denser packings, 
such as thermal equilibration of the system (e.jr via Monte Carlo methods to solve the 



"adaptive shrinking cell" optimization problem 5l|, |52j) or relaxation under pair potentials 
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known to favor high-density configurations. Another possibility would be to combine the 
strengths of the TJ algorithm with those of other lattice packing methods. The ability of 
the TJ algorithm to quickly generate extreme lattices (the inherent structures) could be 
used as a starting point for an algorithm that performs an exhaustive search in the space 



of perfect lattices [20|. Moreover, its efficiency in finding locally-densest lattice packings 
from arbitrary initial conditions could be used to rapidly obtain such packings starting from 
intermediate-density packings generated using other methods [19f . As d increases from one, 
the first dimension in which the densest known packing that is not a Bravais lattice (peri- 
odic packing with a multiple-particle basis) is d — 10, which has a basis of 40. Since the TJ 
algorithm was successfully used to obtain the densest known packings for d < 6 with a large 
multiple-particle basis (up to a basis of 729 for d = 6) [32], it would be interesting to explore 
whether the TJ algorithm could be used to discover currently unknown denser non-lattice 
packings in 10 dimensions or higher. 

For d > 17 dimensions, the TJ algorithm mainly produces lattices that have both a low 
packing density and a minimal kissing number, while still being locally densest, revealing 
a richer and more complex density landscape than in most dimensions less than 17. This 
phenomenon could possibly be exploited to quickly generate low-density extreme lattices in 
very high dimensions. Since these lattices are strictly jammed and have the minimal kissing 
number to ensure mechanical stability, they can be considered to be the lattice analogs to 
the maximally random jammed packings (disordered local-maxima inherent structures) that 



have been generated using the TJ algorithm with many particles per fundamental cell [32 ]. 
Such configurations could be generated in much higher dimensions than those considered in 
this paper, since the requirement of reaching the ground state would be removed, and the 
TJ algorithm is less resource-intensive when generating suboptimal kissing configurations 
(through the reduced number of constraints). 
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Appendix A: Lattice definitions 

In this appendix, we define some common lattices, following the notation and nomencla- 



ture used in Refs. 
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and 
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The hypercubic Z d lattice is defined by 

Z d = {(xi, . . . , x d ) : Xi GZ} for d > 1 (Al) 

where Z is the set of integers (. . . — 3, —2, — 1, 0, 1, 2, 3 . . .) and Xi, . . . , Xd denote the compo- 
nents of a lattice vector. The kissing number of Z d is 2d. A d-dimensional generalization of 
the face-centered-cubic lattice is the checkerboard D d lattice defined by 

D d = {(x u ...,x d )eZ d :x l -\ Vx d even} for d > 2. (A2) 

Its kissing number is 2d(d — 1). Note that D2 is simply the square lattice Z 2 . Another 
generalization of the face-centered-cubic lattice is the root lattice Ad, which is a subset of 
points in Z d+1 , i.e., 

A d = {(x ,xi,...,x d ) G Z d+1 : x + x x H h x d = 0} for d > 1. (A3) 

The kissing number of A d is d(d + 1). In three dimensions, D% and A3 are identical, but D d 
and Ad are inequivalent for d > 4. Another set of root lattices is denoted E d , for d = 6, 
d = 7, and d = 8. The root lattice E 8 is equal to the union of D s and the translation of 
D 8 by (5,|,|,5,|,|,5,5)- The root lattice £7 is the section of E 8 where the sum of the 
lattice coefficients is set equal to zero, and the root lattice Eq is the section of E 7 where the 
sum of the first and eight coefficients is also set equal to zero. Alternatively, vectors in Eg 
perpendicular to any ^-sublattice in E 8 also form Eq. 

The laminated lattice A^ is constructed by stacking layers of a (d — l)-dimensional lam- 
inated lattice A^_! as densely as possible such that the shortest vector in A^ is of equal or 
longer length than the shortest vector in A^_i. This definition does not uniquely define A^ 
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for all dimensions. For d — 11, d — 12, d — 13, and d > 25, there exist multiple laminated 
lattices of equal densities, which we distinguish using superscripts. Many of the previously 
defined lattices are also laminated lattices. For example, A x = Z 1 , A 2 = A 2 , A 3 = D 3 , 
A 4 = D 4 , A 5 = D 5 , A 6 = Eq, A7 = E 7 , and A 8 = E 8 . A particularly interesting laminated 
lattice is the 24-dimensional Leech lattice A 2 4- Finally, the Coxeter-Todd lattice K 12 can be 
defined for 18 dimensions: 

K 12 = {(xu, ••■ , x m , x 2 i, • • • , £ 2 6, x 31 , ■ ■ ■ , x 36 ) : Xij e Z} , (A4) 

where xn~ denotes the components of a lattice vector, subject to the following conditions 

xn + x i2 + x i3 = ie{l,---,6}, (A5) 

xu — Xji = x i2 — Xj 2 = x i3 — Xj 3 mod 3 i, j E {1, • • • , 6}, and (A6) 

xik + x 2k + x 3k + x 4k + x 5k + x 6k = mod 3 k e {1,2,3}. (A7) 

This lattice can be generalized to other dimensions in the range 6 < d < 18 by requiring that 
Kd is the densest section of i*Q+i which either contains or is contained in K\ 2 and taking 
^18 = A 18 . 
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